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Abstract 

In this work wc propose upscaling method for nonlinear Forchheimer flow in 
highly heterogeneous porous media. The generalized Forchheimer law is con- 
sidered for incompressible and slightly-compressible single-phase flows. We use 
recently developed analytical results Jj and write the resulting system in terms 
of a degenerate nonlinear flow equation for the pressure with the nonlinearity 
that depends on the pressure gradient. The coarse scale parameters for the 
steady state problem are determined so that the volumetric average of velocity 
of the flow in the domain on flne scale and on coarse scale are close enough. 
A flow-based coarsening approach is used, where the equivalent permeability 
tensor is flrst evaluated following the streamline of the existing linear cases, and 
successively modified in order to take into account the nonlinear effects. Com- 
pared to previous works [21 13] , our approach relies on recent analytical results 
of Aulisa et al. [T] and combines it with rigorous mathematical upscaling the- 
ory for monotone operators. The developed upscaling algorithm for nonlinear 
steady state problems is effectively used for variety of heterogeneities in the do- 
main of computation. Direct numerical computations for average velocity and 
productivity index justify the usage of the coarse scale parameters obtained for 



'Corresponding author 
Email addresses: eugenio.aulisaOttu.edu (Eugenio Aulisa), 
lidia.bloshaiiskaya@ttu.edu (Lidia Bloshanskaya), efendiev9math.taiiiu.edu (Yalchin 
Efendiev), akif . ibraguimovOttu. edu (Akif Ibragimov) 



Preprint submitted to Elsevier 



March 21, 2013 



the special steady state case in the fully transient problem. Analytical upscaling 
formulas in stratified domain are obtained for the nonlinear case. They correlate 
with high accuracy with numerical results. 

Keywords: upscaling, heterogeneity, Forchhcimer flow, nonlinear flow, 
permeability, productivity index 



1. Introduction 

In recent years, using near well data, e.g., core data, engineers create in- 
creasingly complex and detailed geocellular models which result in highly het- 
erogeneous geological parameters of reservoirs. Such descriptions typically re- 
quire a high number of computational cells which are difficult to solve, e.g., in 
well optimization problems and history matching. To reduce the computational 
complexity, some type of coarsening and upscaling procedures are needed. The 
geological parameters, such as permeability or transmissibility and porosity, 
should be upscaled for each coarse-grid block. 

The variety of approaches for upscaling and multiscale methods of fine scaled 
geological parameters have been proposed for the linear Darcy case (e.g., [31[51 
m [U [HI Un])- These approaches include upscaling methods, e.g., [THl [S] 
and multiscale methods [2l E |5]. In both approaches, a goal is to represent 
the solution on a coarse grid where each coarse-grid block consists of a union 
of connected fine-grid blocks. In upscaling methods, the upscaled permeability 
is calculated in each coarse-grid block by solving local problems with specified 
boundary conditions and calculating the average of the flow flux. Local problems 
can be solved in extended domains for computing the effective properties. In 
multiscale methods, the local multiscale basis functions are computed instead 
of local effective properties and these basis functions are coupled via a global 
formulation. 

The extensions of these methods to nonlinear flows, such as Forchheimer 
flow, are carried out in several papers, see [3 [TT] which are closely related to 
our work. In [3] , the authors consider the use of iterative upscaling techniques 
where at each iteration, local-global upscaling technique is used. In pT], the 
authors use special nonlinear forms for upscaled Forchheimer flows that simplify 
the upscaling calculations. In the current paper, our goal is to carry out rigorous 
nonlinear upscaling using new formulations of Forchheimer flows. 

In current paper, we utilize recent flnding [1], where Forchheimer equation 
is written in an equivalent form using monotone nonlinear permeability func- 
tion depending on gradient of pressure. This equivalent formulation reduces the 
original system of equations for pressure and velocity to one nonlinear parabolic 
or elliptic equation for pressure only. The ellipticity constant of this equa- 
tion degenerates as the pressure gradient converges to inflnity. The rate of the 
degeneration is effectively controlled by the order of Forchheimer polynomial 
and the structure of the coefficients has the important monotonicity proper- 
ties [T], Proposition III. 6 and Lemma III. 10. It allows to prove results on the 
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well-posedness of the initial boundary value problem and apply numerical ho- 
mogenization theory. 

In this paper we present the upscaling algorithm for fluid flow in incom- 
pressible media for two types of fluids, incompressible and slightly compressible. 
Steady state problem for incompressible flow reduces to the degenerate elliptic 
equation, however the corresponding problem for compressible fluid reduces to 
time dependent degenerate parabolic equation. 

In this report we first introduce and investigate the upscaling procedure for 
the time independent problem in case of incompressible fluid. In case of time 
dependent problem the question one should address is that while the solution is 
time dependent, the upscaled parameters are time independent for incompress- 
ible media. We use the upscaled parameters obtained for steady state case in 
the time dependent problem. This procedure is justified by the results obtained 
in our papers [12] and [13] and the numerical experiment presented in this ar- 
ticle. Namely, we will relate the fine scale fully transient solution to the special 
pseudo steady state (PSS) solution. This solution has a form At + W{x), where 
A is a constant and W{x) is a solution of auxiliary steady state boundary value 
problem for the equation with non zero RHS. According to our results in [T^ 
and [13] under some assumptions the pseudo steady state pressure and velocity 
serve as pseudo attractors for fully transient pressure and velocity. To upscale 
the steady state equation we determine the coarse scale porosity and nonlinear 
permeability, so that the average volumetric velocity of the flow is preserved. 
To evaluate the described method for the time dependent case we compare the 
productivity index (PI) of the well on the fine and coarse grids. The PI is 
inversely proportional to the difference between the average of pressure in the 
reservoir and on the well. We select the PI as a criteria for the evaluation of the 
upscaling method as it is widely used by the engineers [131 [TS] . In the numerical 
examples we calculate the difference between the values of the Pis on fine and 
coarse grids. Our numerical results show that the proposed algorithm provides 
accurate results for different heterogeneities and nonlinearities in steady state 
case. Resulting transient velocity and PI on coarse scale also provide accurate 
approximation of corresponding transient parameters on fine scale for heteroge- 
neous fields considered in the paper. We expect the accuracy of the proposed 
method depends on heterogeneities as in a single-phase upscaling {"T^), i.e., for 
highly heterogeneous fields, the accuracy of the method will deteriorate. In this 
paper, our main goal is to propose a method to handle the nonlinearities and, 
thus, we do not consider very highly heterogeneous fields (|16]). 
The paper is organized as follows. In Sec.|2.1|we introduce g-Forchheimer equa- 



tions, review their properties and formulate the problem. In Sec. 2.2 we obtain 



the form of the coarse scale equation for generalized Forchheimer flow. Sec. 2.3 



presents the discussion of convergence results for the transient velocity and PI 
in case of slightly compressible flow and usage of the upscaled parameters from 
the steady state equation in transient case. Sec. [3] is devoted to description of 
upscaling algorithm. In Sec. [4] we obtain the explicit analytical upscaling for 
mulas in case of incompressible fluid for stratified region. In Sec. |5.1| and 5.2 



we 



present the numerical results for the incompressible and slightly compressible 
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flows correspondingly. 



2. Problem statement and Preliminary results 

2.1. Generalized Forchheimer equation 

Darcy equation describes the linear dependence of velocity u on the pressure 
gradient Vp 

u = -lk{x)\/p. (1) 

Here k{x) is symmetric positive definite permeability tensor, /i is the viscosity 
of the fluid. 

Forchheimer equation [17 is known to generalize Darcy's equation to take 
into account inertial terms and has been introduced in the literature in several 
forms. E.g., 

Two term law: u + /3{x)\\u\\u = — ^A;(a;)Vp, 

Three term law; u + ai{x)\\u\\u + a2{x)\\u\\'^u = —j^k{x)\7p, (2) 
Power law: u + bi{x)\\u\\"^^^u ^ ~j^k{x)Vp, 1.6 < m < 2. 

Coefflcients /3{x), ai{x), a2{x) and bi{x) are empirical. 
All these relations can be written in a compact form as 

g{\\u\\,x)u = -l^kix)Wp, (3) 

for some function g{s, x) > for s > 0. We will refer to ([S]) as g -Forchheimer 
{momentum) equation. For simplicity from now on we assume the viscosity 
11 = 1, i.e. j^ = k. 

To develop rigorous numerical homogenization concepts for Forchheimer 
flow, we use the results in [1] which allows writing ([s]) as a monotone rela- 
tion for Vp. Moreover, this allows obtaining the well-posedness results of the 
corresponding initial boundary value problem and allows estimating the residual 
error in numerical homogenization because of monotonicity. It was shown in [T] 
that the monotone relation between velocity and gradient of pressure exists for 
general functions g{s,x) in the form 

k 

g{s,x) = l + ^aj(a:)s"^- = 1 + ai{x)s°'^ + a2{x)s°'^ + . . . + ak{x)s°"' , (4) 

where fc > 0, the exponents satisfy < aj < Ofj+i, and the coefflcients a,j(x) > 
0, j = 1, . . . , fc. Thus deflned function g in (|3| includes all the known cases of 
Forchheimer flow ([ij and ([2|. 

We can define the inverse function 

G{C,x)= J e>0, h{s) = sg{s,x), s>0. (5) 
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We then can obtain the equivalent form of Eq. ([3| 

-G{\\k{x)Vp\\;x)kix)Vp, (6) 
which we caU generalized (nonlinear) Darcy equation. 

Remark 2.1. In the particular case of two-term Forchheimer law, the nonlinear 
permeability tensor G can explicitly be written 

G(^;x)^ (7) 

The g-Forchheimer equation written in the form ^ allows reducing the 
dynamical system to single nonlinear equation of pressure. Namely, we consider 
the continuity equation 

cl^{x)f^=-V-{pu), (8) 

where p is the density of the fluid, and is the rock porosity. For incompressible 
fluid (p = const) , ([s]) reduces to V • u = and combined with the flow equation 
([6| results in the degenerate elliptic equation of pressure only for steady-state 
flow 

V • (G(||fc(a;) Vp||;x)fc(a;) Vp) = 0. (9) 

For slightly compressible fluid (such as the compressible liquid) the equation 
of state takes the form, see il4j , 

p{p) = poe-'P, (10) 

where 7 is the inverse of the compressibility constant. 
Substituting ^ in ^ we get 

0^ = -^V-ii + Wp. (11) 

For slightly compressible fluids 7 is of order 10^^, thus we drop the second 
term in RHS of equation (11 1. Combining it with ([6| we obtain the degenerate 
parabolic equation for pressure 

70(a;)^ =V-(G(||fc(x)Vp||;a;)fc(a;)Vp). (12) 

Equations ^ and ( 12 ) describe the fluid flow on the fine grid. Our aim is to 
devise an upscaling algorithm for the parameters k{x), 4>{x) and Gdj/cVpH ; x) 
in equations ^ and ( 12 ) and obtain the corresponding coarse scale equations. 

2.2. Coarse scale equation in case of incompressible fluid 

To obtain the coarse scale equation for incompressible case, we first rewrite 
Eq. ([9]) in each coarse block fie in a form: 

V -KikVpix) = in n^, (13) 
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where K{ri;x) = GdlryU; a:)??. We assume p — ^ ■ x on dflc, where vector ^ = 
{£,i,£.2) (and ^ — {^1,^2,^,3) in 3D). We then solve Eq. (13l in each coarse block. 
In each coarse block fi^ we define 



(/) = 



Here 

1 

l^cl JO, 

is the volumetric average of the function over fi^- 

We would like to find the upscaled tensor k* and scalar G* , depending on 
fc*, so that 

Then, the upscaled equation takes the form 



V-(G*(Vp*)fc*Vj3*)=0. 



(14) 



It follows that the coarse scale function G* depends on the vector k*'Vp* , 
while the fine scale function G depends on the scalar \\k{x) Vp||. 

2.3. Coarse scale equation in case of slightly compressible fluid 

Unlike the steady state Eq. ^ for incompressible fluid, Eq. ( 12 1 for slightly 
compressible flow is transient in time. Parameters fc, G and on the fine scale 
are, however, time independent. We want to find the upscaled parameters fc*, 
G* and <&* on the coarse scale which are time independent as well. It is difficult 
to use the original equation ( 12 ) for the upscaling procedure directly. Instead 



we will relate the fine scale transient pressure and velocity to the special pseudo 

Teur,. 



steady state solution of Eq. ( 12 ) which will be defined below 



Let U be the domain with the boundary T consisting of two parts F 
The no-flux condition is imposed on Eg 

wi/|r^=0; (15) 

and prescribed total flux condition is imposed on Ti 



vds^Q{t), (16) 



where u is the velocity as in ([6]) and v is the outer normal to the boundary E. 
In [T] it was proved that there exists a special solution pa{x,t) of equation 



(12 1 with boundary condition (15 1 on Ee such that 

dp. 
dt 



const. 



-A 



for all t. 



(17) 



Such solution is called Pseudo Steady State (PSS). Erom definition (17) of PSS 



solution it follows that the corresponding production rate is constant Q{t) 
Q = A\U\ — const, and Ps{x,t) can be written as 



Ps{x,t) 



'7|||i + VE(x) 
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where W{x) is called a basic profile and is a solution of the steady state BVP 

V • iG{\\kix)VW\\;x)kix)VW) = -7 J^^^i^)^ (18) 
us{x) ■ = 0, (19) 
W\r^ = M^), (20) 

with given function (po{x) and constant C. Notice that Vps — WW and the 
corresponding PSS velocity 

u,{x) = -G{\\kVps\\;x)kVp, = -G{\\kVW\\;x) kVW (21) 

is time independent. 

The steady state BVP ( 18 )-([20| will be used to find the upscaling parameters 
for fully transient equation (|12|. On coarse scale the steady state Eq. ( 18 ) will 
take the form 

V-(G*rVp*)^-|^$*. (22) 

The upscaling algorithm for k* , G* and <I>* follows, as previously, a procedure 
of equating average velocities (cf . |5] ) . Consequently, 

II {us)^ — {ul)^ II is sufficiently small. (23) 

Here 



^'u - T77T ; ("s)c/ = T77T '^^^ (24) 



where and it* are the steady state velocities on fine and coarse scales, corre- 
spondingly. 

If initial data is not of basic profile then solution p(x, f) of the original equa- 
tion (12 1 and the corresponding velocity u{x,t) are time dependent. Thus in 
order to justify the upscaling criteria ( 23 1 for general case one should prove con- 
vergence of the corresponding time dependent quantity to the time independent 
one. This property was obtained in |12j and |13j under certain conditions on 
the boundary data. Namely, let 

■ip{x,t)^ p{x, lr. - ||!:| y Pi^^ 

and 

if{x) = ipo{x) - 1^ X 'Poix) ds 

be the deviations from the average on the boundary P^ of the trace of transient 
solution p{x, t) and basic profile W{x) correspondingly. 

We proved that if the differences Q{t) — Q and ip{x, t) — ip{x) converge in 
certain sense to zero at time infinity (see [13|, §3.2), then the PSS velocity Us{x) 
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serves as the steady-state attractor for the fully transient velocity u{x, t) with 
any initial data: 



\\u{x,t) - Us{x)fdx 



as t ^ oo. 



This justifies the usage of criteria (23 1 for the upscaling of coefficients k, G, and 
(f) in fully transient problem. 

To evaluate this method we made a comparison between the productivity 
index on coarse and fine scale using the coefficients k* , G* , and <&* on coarse 
scale. Productivity index is routinely used by engineers in estimation of available 
reserves and optimizing well recovery efficiency (see [HI [151 [18] ) . It is defined 
as follows. Let p{x, t) be the solution of BVP in region U for equation ( 12 1 with 



(PT) is defined as the ratio 



boundary conditions (151 and (16). The Productivity Index/ Diffusive Capacity 



Jit) 



Q{t) 



(25) 



where pjj (t) — pp . (t) is a pressure drawdown; and 



1 



1 



p{x, t) ds. 



In case of PSS flow the productivity index J{t) is time independent and 



J{t) = J. 



PSS 



\u\ Ju 



LW{x)dx 



(26) 



We compare the Jpss corresponding to the solution of the fine scale equation 



( 18 ) with the Jpgg corresponding to the solution of the equation ( 22 ) on coarse 
scale . It is numerically proved that the difference \Jpss ~ Jpss is small (see 



Sec. 5.21 



As it has been already mentioned, in general, the productivity index is time 
dependent. As in the case for velocity, it was proved that if the differences 
Q(t) — Q and tp{x,f) — ip{x) converge in certain sense to zero at time infinity 
(see [H], §3.2), then 



\J{t) - Jpss\ 







as t 



oo. 



Thus the coarse coefficients k* , G* and obtained for the steady state equation 
( [l8l ) , can be used to calculate fully transient productivity index on coarse scale. 

Numerical experiment confirms the theoretical findings above. Figures [l] and 
[2] present the time dependence of velocity and the PI of coarse and fine scales. 
The time dependent values are also compared to the PSS values, which are 
constant in time. As it can be seen from the graphs, in the long term the coarse 
scale time dependent velocity and PI calculated using the upscaled parameters 
from the steady state problem provide good approximation of the corresponding 
fine scale values. 



8 



1.2 




D COS 0.1 0,15 0,2 0,25 

Time 



Figure 1: Time dependence of the average velocity on fine and coarse scales: I 
- PSS (us) on the fine scale; II - PSS (u^) on the coarse scale; III - {u{t)) on the 
coarse scale; IV - {u{t)) on the fine scale 
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Figure 2: Time dependence of the Productivity Index on the fine and coarse 
scales: I - PSS PI on the fine scale; II - PSS PI on the coarse scale; III - PI(i) 
on the fine scale; IV - PI(t) on the coarse scale 



3. Numerical upscaling algorithm 



In this section we present the numerical upscaling algorithm for the steady 
state equations ( 18 ) and We consider 2D rectangular region Q, with horizon- 
tal size Li and vertical size L2 and two orthogonal grids: fine x — (xi, a;2)-scale 
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and coarse X = (Xi, X2)-scale (see Fig. [s]). 

^^^2 




Figure 3: Fine and Coarse Scale 



The porous media on the fine scale is considered to be isotropic, and perme- 
ability tensor k{x) is a scalar function k{xi,X2)- The fine scaled equation ( [l8| 
or ([9]) with parameters k{xi^X2), 4>{xi,X2) and G(||fcVp|| ; xi, 2:2) is upscaled to 



the coarse scale equation (22) or (14 1 with parameters k*{Xi^X2), $*(Xi,X2) 
and G*{Vp*;Xi,X2) so that condition ([23| is satisfied. The Us and u* in ([23| 



are the velocity on the fine and coarse scale correspondingly. 

Our approach is purely local, so the algorithm is described for single coarse 
block with boundary dVic- For simplicity we take Vic to be the rectangular 
[0, h] X [0, h] with the area \Vlc\ — h ■ h, see Fig. [sj 

For each coarse block the two-step procedure is performed: 



Step 1: the equivalent permeability tensor k* 
methods; 



is obtained using linear upscaling 



Step 2: the equivalent nonlinear coefficient G* is obtained using k* . 

Step 1. Procedure to obtain k* . In order to obtain full permeability tensor 
k* we use the standard local procedure via volume averages of velocity and 
pressure gradients, see for example [Ul. We solve two fiow problems in each 
block with periodic boundary conditions. Namely, let pi and p2 be the solutions 
of the fine scale equation in coarse block fie 



V • {k{xi,X2)Vp) = 

with boundary conditions: 
' pi{xi,0) = pi{xi,l2) for a;ie[0,/i] 



(27) 



ui(xi,0) • = 
Pi(0,a;2) = 0; 



-Ui{xi,l2) ■ Vi] 

Pi{h,X2) = 1; 



P2iO,X2) ^ P2{h,X2) for a;2e[0,/2]; 
"2(0, 2:2) • vi = -1*2(^1,2:2) • 1^2; 
P2{xi,0)^0; ^2(2:1,^2) = 1- 

(28) 
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Here u,; is the velocity vector corresponding to the pressure distribution pi, 
i = 1,2. 

The four elements of the upscaled permeability k* are then calculated from 
two vector equations: 

(u,) = -k* (Vp,) , 1^1,2. 

The upscaled porosity <i>* is computed via integral averaging on the coarse 
block following classical approach, e.g. |19) : 

$* - (0(xi,X2)) = (f>{x,,X2)dn,. (29) 



Step 2 : Procedure to obtain G* . We use the upscaled permeability k* 
to determine the nonlinear coefficient G* via pure local averaging. As it was 
mentioned in Sec. 2.2 unlike the fine-scale function G depending on ||fcVp||, the 
upscaled G* depends on the vector Vp* itself. Let f = {£,1,^2) be the gradient 
of pressure in coarse block Vic- For fixed ^, the G* is a constant. If p* is the 



solution of coarse scale Eq. ( 14 ) with boundary condition 



P Ian, 



£.1X1 + £2X2, 



(30) 



then it is also the solution of equation V • (fc*Vp*) — with the same boundary 
condition. 

We determine G* so that 

II II -IK") II, (31) 

where u* = —G*k*Vp* is the velocity on coarse scale and u = —GkVp is the 
velocity on fine scale corresponding to the solution p{x, t) of ^ with boundary 
condition p|ao, = ^1X1 + £,2X2- Then for fixed Ci,'?2 we have: 



G*(6,e2) 



\k* (Vp* 



(32) 



Using formula ( 32 1 we numerically construct the table of values of G* for 
6,6 e (-00,00). It follows that G*(0,0) = 1, G*(6,6) 1 if 1161 ^ 0, 
G*(6-6) — > if 1161 00 and G* possesses certain symmetry: G*(6:6) = 
G*(— 6,— 6)- It is thus sufficient to consider 6 S (—00,00) and £,2 > only. 
It is worth mentioning that the special attention should be paid to the way 
the domain for the ^ is discretized. Taking the grid to be too fine makes the 
calculations overly expensive, however the sparse grid does not allow to capture 
the features of nonlinearity of the process. We will use the non uniform grid, 
where the subsequent point is calculated on the basis of the deviation between 
the preceding values of the function. Namely, 6 ^'^.d 6 taken from the set 
7y„, n = 0, 1, 2, . . . , 
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where the first three values are taken a priori: rjQ — 0, and rji, 772 to be small 
enough. Next value jyn+i is chosen so that 



71 + 1 



< £„, where = G* {r]m,rim), m ^ n,n + I 



for some set value e„. The stopping criteria for the computation is G* < e and 
|(G* )'| < Ed so that the value of G* as well as its variation are sufficiently small. 
The shape of function G* is presented on Fig|4j 




-200 



Figure 4: The shape of function G*(^i,^2) 



4. Analytical Upscaling for the Layered Porous Media 

Here we will present the analytical upscaling formula for the nonlinear Forch- 
heimer flow of incompressible fluid in layered porous media. Consider a rect- 
angular region R of horizontal size L and vertical size H. The region has a 
horizontal multilayer structure, and is composed by n layers, see Fig. [5] Each 
ith layer, i = 1, . . . , n, has vertical size hi and is characterized by constant iso- 
topic permeability ki and f;-Forchheimer polynomial g{s,x) — gi{s) with con- 
stant coefficients or, equivalently, by the nonlinear function Gi — Gi{ki\\Vpi\\). 
We assume that the type of nonlinearity is the same for each layer, while the 
coefficients of ^-polynomials can be different. 

Under these assumptions, in each ith layer, i = 1, . . . , n, equations ([3|, (|6| 
and ^ yield 

gi{\\ui\\)ui = -kiVpi, (33) 
Ui = -Gik^Vpi, (34) 
V • (G,fc,Vp,) = 0. (35) 

Here Ui, pi and Qi are, correspondingly, velocity, pressure and the total bound- 
ary flux in ith layer Ri and || • || is the I2 vector norm. 
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^2, G2 



h2 



kni Gn hr, 
I 

Figure 5: Layered porous media, region R 



We assume that flow within the whole block R is subject to the equation with 
the same type of nonlinearity as in each layer. We aim to find the equivalent 
homogeneous block permeability k* and nonlinear coefficient G* — G*{k*Vp*) 



for two types of flow: flow parallel to the layers (Sec. 4.1 ) and flows perpendicu- 



lar to the layers (Sec. 4.2 ). The upscaled parameters are determined so that the 
total flux of the system stays the same as with nonhomogeneous parameters. 
The comparison between the obtained analytical results and numerical compu- 



tations using the method in Sec.[3]are presented in Sec. 5.1 Note, that the case 



when Qi = const, is the same as Darcy case and the upscaling formulas for k* 
are the same as in [20] . 

Flow Parallel to the Layers 
We impose the following boundary conditions on boundaries of R 

• Pi — P* = Po, on the left boundary, i ~ 1, . . . ,n, 

• Pi — P* = Pi oil the right boundary, i = 1, . . . , n, 

• M • I' = on the bottom and top boundaries, 

where v is the outward normal, and pi > po- Under these conditions the flow 



is parallel to the layers and the solution of Eq. (35 1 is linear in x. The pressure 
gradient is constant and is equal to Vp — 0), where ^ = {pi ~po)/L. In each 
layer the vertical velocity component is identically zero, while the horizontal 



component is constant in each layer and, according to ( |34[ ), is equal 

Ui^ = Ui ^ ~Gi{ki^) ki ^. (36) 

On the other hand the outgoing flux is equal to incoming fiux and is equal to 
the sum of fluxes in the ith layer: 

n n n 

Q = = (37) 

2—1 2—1 2—1 
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The flux is zero on the top and bottom boundaries. 

We now consider the analogous block with the same boundary conditions 
and permeability k* and nonlinear function G* resulting in the same flux Q. In 
this case the flux is 

Q = -u*H = G*k*H^, 

where u* is constant horizontal component of upscaled velocity. Expression for 
Q above and (37) yield 



71 1 ^ 

Gjkihj ^ G*k*H or u* ~ — Ujh. 



(38) 



First we consider the limiting linear Darcy case Gi — 1. In this case G* — 1. 
We then find an expression for k* 



1 " 



i=l 



In view of ( 39 ) the general expression for G* follows from ( 38 ) 



G* 



1 1 



^ ^ Giklhi 



i=l 



X^j— 1 kihi 



(39) 



(40) 



Formulas ( 39 1 and ( 40 1 can be generalized in case when the parameters k 
k{x2), G — G{k(^,X2) are continuous functions: 



k* = 



H 



H 



k{x2) dx2; 



G* = 



1 1 



H 



G{k^,X2)k{x2) dx2 



Alternatively, using g-Forchhcimer equation (33) with \\ui\\ = Ui, the upscal- 
ing formula for the (^-polynomial can be obtained: 



g*{u*) 



From here we can obtain the upscaled coefficients a*, j = 1, . . . ,r7i corre- 
sponding to the power , for the g-polynomial in domain R. In particular in 
case of two-terms law as in Remark 2.1 the upscaled Forchheimer coefhcient can 
be obtained explicitly in the form 



,2 hi 



H 



(41) 



where Pi is coefficient corresponding to ith layer. 

The coefficient /3* depends explicitly on ^. The two limiting cases are 



lim B* 



14 



and 



lim / 



yi^i=l y H J \l^i=l y 13, H 



In case when the parameters k — k{x2) and j3 — f3(x2) are continuous 
functions, the expression (41) for /3* yields 



1 rf/3(a;2) ^ 



2 



-^0 l + Vl+4/3(2;2)fc(a;2)? 



Flow Perpendicular to the Layers 

Let consider the same geometry and let impose the following boundary con- 
ditions 

• Pi\x2=o — P*\x2=o = Pa, on the bottom boundary, 

• Pi\x2=H = P*\x2=H = Pn, on the top boundary, 

• u • I' = on the left and right boundaries. 

In this case the flow is perpendicular to the layers and the horizontal velocity 
component Ui-^ is identically zero, while the vertical component of velocity = 
u is constant in each layer. The pressure gradient in ith layer is equal to 

Vpj==(0,^,), where = 

where pi is the pressure measured at the top of ith layer for i = 1, . . . ,n. Then, 
according to (34), vertical component of velocity is equal to u = —Giki^i. 

It thus follows that the flux is constant and in each layer is equal to Q = 
—uL = Giki^iL. Then the pressure gradient in ith layer is 

e. = ^v^, z = l,...,n. (42) 

We again want to identify the equivalent homogeneous parameters k* and 
G* in the region R resulting in the same flux Q* = Q. The pressure gradient in 
the domain R is 

?-^4EeA. (43) 

We get the expression for the flux 

Q = -u*L = G*k*^L. (44) 
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Plugging ( 43 ) in ( 44 ) and using ( 42 ) we get 

1 

4=1 



G*k* 



1 " 

-Y 



and it follows: 



1 1 " 
k* H ^ 



K 
k, 



1 

G* 



k ^ — ^ 1 
H^^Giki 



(45) 



Alternatively, using g-Forchheimer equation (33) gi{u)u = —hS^i, the up- 
scaling formula for the ^-polynomial can be obtained: 

, k*^ , M 



i=l fc. 



and thus 



k ^ T /i^ 

i7 Z.«^-'^ 



J=l,-. 



(46) 



where aj^i is the coefficient of g-polynomial corresponding to power s"^ , j = 
1, . . . , TO (see Eq. Q) in ith layer. 

In case when the parameters k — k(x2), G ~ G(ki^(xi)^xi) and aj = 0^(2:2), 
J = 1, . . . ,TO, are continuous functions in x^-, Eqs. (45) and (|46|) yield 



1 



1 



1 



fc* 



dX2 



i^7o K^^y G* H Jo G{ki{x2),x^)k(x2y ' 



k*_ aj{x2) 
H Jo k{x2) 



dx-2. 



5. Numerical Results 



In this section we numerically illustrate the upscaling algorithm described 
in Sec. [3] for the incompressible and slightly compressible fluids. The considered 
cases of the permeability distribution on fine scale are presented in Fig. [6] The 
obtained upscaling errors are relatively small, since we did not consider large 
heterogeneities, but instead focused on the upscaling method for the nonlinear 
flow. 



5.1. Numerical Results for Incompressible Fluid 

In this section we present the numerical results for incompressible flow. Sev- 
eral approaches are compared: the upscaling algorithm Sec. [3] and the analytical 
formulas obtained in Sec. ID 

On the fine scale the pressure is subject to equation ^ in the region — 
[0, Li] X [0,^2] with the boundary conditions 

p(0,X2)=0; p{L^,X2)^l; ||(xi,0) = ^(a;i,L2)=0. 
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We report the relative error in the averaged velocities originated from the 
upscaling of equation ^ to (14): 

II (u) - (u*) II 



II W II 



(47) 



The errors for the layered system are reported in Table [T] (flow parallel to the 
layers, permeability as in Fig. 6b) and Table [2] (flow perpendicular to the layers, 
permeability as in Fig. 6a). The results for the case of random system are 
reported in Table [3j 

For each case we compare results obtained in three different ways: 1) an- 



alytical formulas Eqs. (39), (40) for flow parallel to the system (denoted by 
"Av 



in the tables); 2) analytical formulas Eqs. ( |45| for flow perpendicular 
to the system (denoted by "Av tl"" in the tables); 3) numerical approach de- 
scribed in Sec. |3] (denoted by "Num" in the tables). For the layered system, the 
corresponding analytical formula gives the exact result. 

The calculations are performed for different orders of nonlinearity: 

1. linear Darcy case, (3 = 0; 

2. two-term Forchheimer law, with nonlinear function G as in (0. In this 
case the coefficient (3{x) in ([7| is taken with its relative magnitude Al3/(3min = 
1, 10, 100, where A/3 is the difference between the maximum value of f3{x) and 
the minimum value f3min- 

The coarse grid is considered to be 20 x 20 where each of the coarse-grid 
block contains 20 x 20 fine blocks, relative magnitude of the permeability is 
II Afc||//c,„i„ = 10, where Ak is the difference between the maximum value of 
k{x) and the minimum value kj^in- For the layered system we consider both 
fine and coarse grids to be square. For the random system three different cases 
are considered: ^ — 0.1,1,10. Here Hi,H2 and /ii,ft.2 are the size of 

coarse and fine cells correspondingly. 

Both analytic averaging formulas are computationally cheap. They show 



different performances: formulas (39), (40), derived for flow parallel to layers 



are consistently better than formula ( 45 ) derived for flow perpendicular to layers. 



From Table [2] it can be seen that the accuracy of both formulas decreases as the 
relative magnitude of nonlinear coefficient 



increases. 





(3 = 




^13 


-100 


Av ^ 














Av fr 





5.36e-4 


2.36e-3 


1.7e-2 


Num 















Table 1: Upscaling errors for numerical and analytical methods, permeability 
of Fig. [6b 



5.2. Numerical Results for Slightly Compressible Fluid 

Numerical results for upscaling in case of slightly-compressible flow are pre- 
sented in Tables [5][6] The coarse grid is taken to be 4 x 4 where each of them 
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/3 = 






-100 


Av 





8.38e-4 


9.56e-4 


7.6e-4 


Av t|- 














Nuni 















Table 2: Upscaling errors for numerical and analytical methods, permeability 
of Fig. 



6a 



Hi _ hi _ 1 ri 

H2 ~ h2 ~ " 


/3 = 


A/i _ -1 


^'^ = 10 


-100 


Av ^ 


8.5e-3 


7.39e-3 


8.63e-3 


2.43e-2 


Av tr 


8.5e-3 


1.87e-2 


8.44e-2 


0.18e-l 


Num 


8.5e-3 


6.52e-3 


9.65e-3 


1.70e-2 


Hi _ /ii _ 1 


/3 = 


A^i _ 1 


=10 


-100 


Av ^ 


1.38e-2 


1.52e-2 


2.76e-2 


9.81e-2 


Av -f|- 


1.38e-2 


1.35e-2 


6.55e-2 


1.20e-l 


Num 


1.38e-2 


9.53e-3 


1.25e-2 


1.60e-2 


Hi _ hi _ n ^ 
H2 ~ 'l2 — "-^ 


/3 = 


A/; _ -, 


= 10 

;3„„„ 


-100 


Av ^ 


1.86e-2 


2.21e-2 


3.87e-2 


1.32e-l 


Av tr 


1.86e-2 


l.lle-2 


5.59e-2 


9.33e-2 


Num 


1.86e-2 


1.09e-3 


1.54e-2 


1.99e-2 



Table 3: Upscaling errors for numerical and analytical methods, randomly dis- 
tributed permeability 



consists of 64 x 64 fine blocks. 

On the fine scale the pressure is subject to equation ( 18 ) in the region 17 = 
[0, Li] X [0; L2] with the boundary conditions The following boundary conditions 
are imposed 

p{L^,X2) = 0; ^^(O,^^) = J^(a:i,0) = ^^(xi,^^) = 0. 
We report the relative errors in the average velocity and the PI between exact 



and upscaled solution, given by (47) and \PI — PI*\/PI. 

The calculations are performed for linear Darcy case and two-term Forch- 
heimer law, with G{xi, X2) as inj7|. Four distributions of the fine permeability 
field k are considered, see Fig. I6j The commonly used empirical formulas to 
relate porosity and permeability k are of the form 

where ao = 0.2222 (Wilhe and Rose equation), ag = 0.2272 (Timur equation), 
ao = 1/3 (Morris and Biggs equation), ao = 0.1961 (Berg equation). In our 
calculations we take 4> = 0.1 • fc" where a — 0.33,0.25,0.2. Following [5T] the 
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Forchheimer coefficient (3 is taken 




Numerical results show that the proposed upscaling algorithm provides small 
errors for the upscaled average velocity and productivity index. In particular, 
the relative errors are less than 5 % in all cases. Though, we use a different 
expression for the case of compressible flow compared to incompressible flow, 
we observe that the velocity errors become larger for nonlinear flows. 





a = 1/3 


a = 1/4 


a = 1/5 




PI err 


Vel err 


PI err 


Vel err 


PI err 


Vel Err 


Darcy 


2.93C-2 


5.5C-3 


2.86C-2 


4.3C-3 


2.81e-2 


3.5e-3 


2-Forch 


2.4e-3 


7.6e-3 


3.5e-3 


7.1e-3 


4.3e-3 


7.2e-3 


Table 4: Upscaling errors, permeability of Fig. 


6c 


I 




a = 1/3 


a = 1/4 


a = 1/5 




PI err 


Vel err 


PI err 


Vel err 


PI err 


Vel Err 


Darcy 


1.3e-3 


6.4e-3 


1.5e-3 


7.4e-3 


1.6e-3 


8.0e-3 


2-Forch 


3.9e-3 


8.3e-3 


4.4e-3 


9.2e-3 


4.7C-3 


9.9e-3 


Table 5: Upscaling errors, permeability of Fig. 


61; 


) 




a = 1/3 


a = 1/4 


a = 1/5 




PI err 


Vel err 


PI err 


Vel err 


PI err 


Vel Err 


Darcy 


3.9e-3 


3.07e-2 


3.7e-3 


3.12e-2 


3.5e-3 


3.15e-2 


2-Forch 


1.72-2 


3.35-2 


1.84e-2 


3.48e-2 


1.92e-2 


3.57e-2 



Table 6: Upscaling errors, permeability of Fig. 
6. Conclusions 



6c 



The developed upscaling algorithm for nonlinear steady state problems 
can can be effectively used for p-Laplacian type equations of the form (|9| 



and ( 18 ) and for variety of heterogeneities in the domain of computation. 



The coarse scale parameters k* , G* and $* are determined so that the 
volumetric average of velocity of the flow in the domain on flne scale and 
on coarse scale are close enough. 

The numerical results show that the proposed method can be used to ap- 
proximate the Productivity Index (PI) of the well in the bounded domain 
on the coarse scale. 

Analytical upscaling formulas in stratified domain are obtained for the 
nonlinear case. 
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a = 


1/3 


a — 


1/4 


a = 


1/5 




PI err 


Vel err 


PI err 


Vel err 


PI err 


Vel Err 


Darcy 


1.98e-2 


5.7e-3 


1.9e-2 


5.5e-3 


1.84e-2 


5.5e-3 


2-Forch 


4.4e-3 


9.7C-3 


5.3e-3 


1.05e-2 


5.9e-3 


l.lOe-2 



Table 7: Upscaling errors, permeability of Fig. 6d 



• In our results for the nonlinear problems, the upscaled parameters depend 
on the range of boundary data. 

• Our results on asymptotic behavior of fully transient velocity and PI and 
actual numerical computations justify the usage of the coarse scale pa- 
rameters fc*, G* and $* obtained for the steady state case in the fully 



transient problem (12). 



A. Appendix 

In the Appendix we present the prospective improvement of the upscaling 
method, which in some cases can give better results with only small increase in 
computational costs (see Tables 9p0). 

Here we consider the upscaled fc* and $* to be of a form 



K21 K22 ■ K 



where K = 



Kq + Ki{xi - Ci) + K2ix2 - C2); 
- Ci) + ^C{X2 - C2) + "^dUxi - Ci f - {X2 - C2) 



(Al) 

(A.2) 
(A.3) 



Here c = (ci, C2) is the central point of the coarse cell Sic and Kii,Ki2, K21, K22, 
Kq,Ki,K2, ^a,^b,^c and are constants to be determined. 

First, the permeability tensor k* is determined, and then it is used to upscale 
the porosity 4>. 



I. Permeability fc*. We first obtain the polynomial K in (A.2 1 as the least 
square approximation of permeability k{xi,X2)- 



Ko 



1 



k{xi, X2) dflc 



K, 



Jj-^^ k{xi,X2){xi - Ci) dQc 



1,2. 



(A.4) 

We now will use the constant elements Kn, K12, K21 and ^"22 of matrix 
fc* to "correct" K{xi,X2) so that the average velocities on the fine and coarse 
scales are the same. For this purpose we will modify our approach presented in 
Sec.m 

We consider two linear fine-scale equations with zero RHS, the exact equation 



(27) and the averaged equation 

V-(Z(xi,X2)VP)=0. 



(A.5) 
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Each equation we solve twice. Namely, let pi and Pi be the solutions of (271 
and (A. 5), correspondingly, subject to boundary conditions (28l) and let p2 
and P2 be the solutions of (27) and (A.5I, correspondingly, subject to boundary 
conditions (28)). 

We now equate the velocity averages on fine and coarse scale, with the coarse 
scale velocity —KVPi, i — 1,2, "corrected" with the elements of the matrix k*: 



1 



1 



i = 1,2. 



(A.6) 



Solving the linear system of four equations (A.6) gives the values of Ku, 

K21, K22- 

II. Porosity $*. To find the coefficients ^a, ^b, and in ( |A.3[ ) we 
consider two equations: the upscaled linear equation with function ( |A.3 ) in RHS 

-V-(fc*VP) = ^A + '^Bixi~Ci) + <i>cix2-C2) + <i>D[ixi-Ci)^-ix2-C2)% (A.7) 

and the exact linear equation 

■{k{xi,X2)Vpp) = (b{xi,X2). (A.8) 
Zero Dirichlet boundary conditions are imposed: 

Pldn^ = pp\dn, = 0. 

Coefficients ^a, ^b, and $d are determined so that the boundary fluxes 
through the faces of the coarse cell corresponding to the upscaled solution P 
and fine scale solution pp are equal. 

Due to the linearity of equation (A.7) solution P is the linear combination 

P^<S>a-Pa + <^b-Pb + ^c-Pc + <^d-Pd, (A.9) 
where P\, X = A,B, C, D, are the solutions of BVPs 

-V ■ik*VPx)^ fx{xi,X2), (A.IO) 
PaIso, =0 

with the RHS 



f\ixi,X2) 



ci 

C2 



for 


A 




for 


A 


= B, 


for 


A 


= c, 


for 


A 


= D. 



{Xi - Ci)2 - {X2 - C2f 

The boundary fluxes through four faces of the coarse cell are then related by 
the same expression as (A.9). Solving the resulting system of four equations for 
4>A, A = A, P, C, D, we obtain the expression for $*. 

The corresponding numerical results are presented in Tables [9p0} For the 
reader's convenience the results presented in Tables [5][6] are included here once 
again for the comparison. We compare errors for the upscaling algorithm with 
A;* and $* calculated using different approaches described above. Four cases are 
considered for each of linear and nonlinear case. 
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(i) k* and $* are calculated as in Sec. [s] (denoted as k* - C; - C); 

(ii) k* is as in ( |A.1|), while <i>* is a constant ([29]) (denoted as k* - P; $* - C); 

(iii) k* is as in Sec. 3] while $* is as in (|A.3| (denoted as k* - C; $* - P); 

(iv) k* is as in (|A.l ) and $* is as in (|A.3[) (denoted as k* - P; $* - P). 



Approach (iv) provides consistently better results for both error in PI and 
velocity. Though it is computationally more expensive, the increase is negligible 
in nonlinear case where the main computational expense comes from obtaining 
the function G*. Approach (ii) is routinely comparable to (iv), and for some 
permeabihty distributions {k as in Fig. 6a a = 1/4, 1/5) is even better. With 
that it only amounts to computation of coefficients ( A.4 ) . As expected, approach 
(iii) is routinely the worst, as there is not sufficient information for correct 
estimation of coefhcients of polynomial (A.3). 
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a = 


1/3 


a — 


1/4 


a = 


1/5 




PI err 


Vel err 


PI err 


Vel err 


PI err 


Vel Err 


Lin. 


k* - 


P; $* - P 


5.5e-3 


L9e-3 


5.6e-3 


2.1e-3 


5.6e-3 


2.2e-3 


Lin. 


k* - 


P;' $* - C 


9.3e-3 


5.5e-3 


8.6e-3 


4.3e-3 


8.2e-3 


3.5e-3 


Lin. 


k* - 


C; $* - P 


2.17e-2 


7.9e-3 


2.17e-2 


8.1e-3 


2.18e-2 


8.1e-3 


Lin. 


k* - 


C; - C 


2.93e-2 


5.5e-3 


2.86C-2 


4.3e-3 


2.81e-2 


3.5e-3 


k* 


- P 


$* _ p 


L6e-3 


6.7e-3 


2.5e-3 


7.2e-3 


3.3e-3 


7.6e-3 


k* 


- P 


$* - C 


1.4e-3 


7.2e-3 


2.3e-3 


6.6e-3 


3.0e-3 


6.5e-3 


k* 


- C 


$* _ p 


2.4e-3 


1.08e-2 


3.5e-3 


1.13e-2 


4.4e-3 


1.17e-2 


k* 


- C 


$* - c 


2.4e-3 


7.6e-3 


3.5e-3 


7.1e-3 


4.3e-3 


7.2e-3 



Table 8: Comparison of upscaling errors, permeability of Fig. 6a 





a = 


:l/3 


a - 


1/4 


a = 


= 1/5 




PI err 


Vel err 


PI err 


Vel err 


PI err 


Vel Err 


Lin. 


k* - 


P; - 


p 


1.8e-3 


2.1e-3 


1.8e-3 


2.4e-3 


1.7e-3 


2.5e-3 


Lin. 


k* - 


P; <i>* - 


c 


2.5e-4 


3.9e-3 


1.4e-4 


2.9e-3 


6.6e-5 


2.3e-3 


Lin. 


k* - 


c'; $* - 


p 


3.6e-3 


1.68e-2 


3.6e-3 


1.70e-2 


3.7e-3 


1.72e-2 


Lin. 


k* - 


C; - 


c 


1.3e-3 


6.4e-3 


1.5e-3 


7.4e-3 


1.6e-3 


8.0e-3 


k* 


- P 


$* _ p 




3.5e-3 


6.8e-3 


4.0e-3 


7.2e-3 


4.4e-3 


7.4e-3 


k* 


- P 


$* - C 




3.7e-3 


7.3e-3 


4.2e-3 


7.1e-3 


4.5e-3 


7.0e-3 


k* 


- C 


$* - p 




3.8e-3 


1.76e-2 


4.3e-3 


1.80e-2 


4.6e-3 


1.82e-2 


k* 


- C 


$* - c 




3.9e-3 


8.3e-3 


4.4e-3 


9.2e-3 


4.7e-3 


9.9e-3 



Table 9: Comparison of upscaling errors, permeability of Fig. 



6b 





a = 


1/3 


a = 


1/4 


a = 


1/5 




PI err 


Vel err 


PI err 


Vel err 


PI err 


Vel Err 


Lin. 


k* - 


C; $* - C 


3.9e-3 


3.07e-2 


3.7e-3 


3.12e-2 


3.5e-3 


3.15C-2 


Lin. 


k* - 


P; $* - C 


5.9-3 


1.27e-2 


6.1-3 


1.23e-2 


6.1-3 


1.20C-2 


Lin. 


k* - 


C; $* - P 


8.4-4 


3.88e-2 


7.7-4 


3.89e-2 


7.2-4 


3.90e-2 


Lin. 


k* - 


P; $* - P 


5.9-3 


1.23C-2 


5.9-3 


1.23e-2 


6.0-3 


1.23C-2 




- C 


$* - c 


1.72-2 


3.35-2 


1.84C-2 


3.48e-2 


1.92e-2 


3.57e-2 


k* 


- P 


$* - C 


l.lOe-2 


1.36C-2 


1.26e-2 


1.36e-2 


1.38e-2 


1.4e-2 


k* 


- C 


$* _ p 


9.6e-3 


3.9e-2 


l.lle-2 


3.9e-2 


1.21e-2 


3.9e-2 


k* 


- P 


$* _ p 


1.49e-4 


2.0e-2 


1.63e-2 


2.12e-2 


1.39e-2 


1.43e-2 



Table 10: Comparison of upscaling errors, permeability of Fig. 6c 
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a = 


1/3 


a = 


1/4 


a — 


1/5 




PI err 


Vel err 


PI err 


Vel err 


PI err 


Vel Err 


Lin. 


k* - 


P; $* - P 


6.2e-5 


4.5e-5 


6.7C-5 


4.7e-5 


7.2e-5 


4.8e-5 


Lin. 


k* - 


P; $* - C 


3.1e-3 


4.5e-3 


2.3e-3 


3.4e-3 


1.8e-3 


2.7e-3 


Lin. 


k* - 


C; $* - P 


1.21e-3 


9.4e-3 


1.21e-2 


9.5e-3 


1.2e-2 


9.6e-3 


Lin. 


k* - 


C; $* - C 


1.98e-2 


5.7e-3 


1.9e-2 


5.5e-3 


1.84e-2 


5.5e-3 


k* 


- P 


$* _ p 


3.5e-3 


5.9e-3 


4.3e-3 


6.4e-3 


4.9e-3 


6.8e-3 


k* 


- P 


$* - C 


3.7e-3 


6.9e-3 


4.5e-3 


6.6e-3 


5.0e-3 


6.7e-3 


k* 


- C 


$* _ p 


4.0e-3 


1.32e-2 


4.9e-3 


1.38e-2 


5.6e-3 


1.42e-2 


k* 


- C 


$* _ c 


4.4e-3 


9.7e-3 


5.3e-3 


1.05e-2 


5.9e-3 


l.lOe-2 



Table 11: Comparison of upscaling errors, permeability of Fig. 6d 
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